Fast and accurate computation of the Fourier transform of an image

We use the Battle.Lemarié scaling function in an algorithm for fast computation of the Fourier transform of a piecewise smooth function f. Namely, we compute for —N m, n N (with a given accuracy ) the integrals * ,1 1 • f(m, n) = I I 1(2;,y) dy (0.1) Jo JO in O(ND)+O(N2 log N) operations, where ND is the number of subdomains where the function f is smooth. We consider an application of this algorithm to image processing. Notwithstanding that it might be advantageous to consider an image as a piecewise smooth function f, it is a common practice in image processing to simply take the FFT of the pixel values of the image in order to evaluate the Fourier transform. Due to the jump discontinuities of the function f, the accuracy of such a computation is poor. We propose our algorithm as a tool for the accurate computation of the Fourier transform of an image since the direct evaluation of (0.1) is very costly. 1. IMAGE AS A PIECEWISE SMOOTH FUNCTION It is natural and useful to consider an image as a piecewise smooth function. For example, the goal of segmentation algorithms is to find the boundaries of smooth subdomains of an image. Furthermore, at the level of pixels, one may consider an image as a collection of tiny squares with different values so that the total image is a linear combination of characteristic functions of elementary squares as in the example in Figure 3. Yet, computing the Fourier transform of a piecewise smooth function is not an entirely trivial matter, especially if the number of discontinuities of f is large or the subdomains are complicated. If we use the FFT of the pixel values of an N x N image (which is equivalent to using the trapezoidal rule in (0.1)) then, in some directions, the error will decay only as 1/N. In other words, instead of a piecewise smooth function, we work with an oscillatory function as in the example in Figure 1. The cost of the direct evaluation of (0.1) for images is prohibitive. The direct algorithm for evaluating the Fourier transform of a linear combination of characteristic functions of elementary squares (or rectangles) computes f(m, n) = E fj(m, n), for —N m, n N, as a sum of contributions from each rectangle [al, b,] x [ci, d1], e2mb1 — e21m(1 e2'1' — e2T'1 fz(m, n) = zi ( —2irim ) ( —2irin ) (1.2) where zi are constants, 1 = 1,.. . , ND and ND is the number of rectangles. Though accurate, such evaluation of the Fourier transform requires O(N2 .ND) operations and since typically for images ND N2, the direct approach is not practical. Thus, there is a need for a fast and accurate algorithm to evaluate (0.1). 244 ISPIE Vol. 2277 0819446010/94/$6.OO